#ifndef INTERPOLATION_HPP
#define INTERPOLATION_HPP

#include <stdexcept>
#include <cmath>
#include <vector>
#include <iostream>

class Interpolation {
public:
    virtual std::vector<double> interpolate(double x) = 0; // 纯虚函数，用于在特定点进行插值

    virtual ~Interpolation() = default;
};

class NewtonInterpolation : public Interpolation {
private:
    std::vector<double> xVec; // 存储插值点的横坐标
    std::vector<double> yVec; // 存储插值点的纵坐标

public:
    NewtonInterpolation(const std::vector<double>& xVec, const std::vector<double>& yVec)
        : xVec(xVec), yVec(yVec) {
        if (xVec.size() != yVec.size()) {
            throw std::runtime_error("Error: The size of xVec and yVec must be equal.");
        }
    }

    std::vector<double> interpolate(double x) override {
        std::vector<double> dividedDifferences(yVec); // 初始化差商表

        // 构建差商表
        for (size_t j = 1; j < yVec.size(); ++j) {
            for (size_t i = yVec.size() - 1; i >= j; --i) {
                dividedDifferences[i] = (dividedDifferences[i] - dividedDifferences[i - 1]) / (xVec[i] - xVec[i - j]);
            }
        }

        // 构建牛顿插值多项式
        double result = 0.0;
        for (size_t i = 0; i < dividedDifferences.size(); ++i) {
            double term = dividedDifferences[i];
            for (size_t j = 0; j < i; ++j) {
                term *= (x - xVec[j]);
            }
            result += term;
        }
        return {result};
    }
};

class HermiteInterpolation : public Interpolation {
private:
    std::vector<std::vector<double>> M;
    std::vector<double> x;
    size_t c,r;
public:
    HermiteInterpolation(const std::vector<std::vector<double>>& M) : M(M) {
        c = M[0].size();
        r = M.size() - 1;
        x.resize(c*r);
        for (size_t i = 0; i < c; ++i) {
            for (size_t j = 0; j < r; ++j){
            x[i*r + j] = M[0][i];
            }
        }

    }

    std::vector<double> interpolate(double x1) override {
        std::vector<double> dividedDifferences(c*r, 0.0);
        int t = 1;
        // 构建差商表
        for (size_t i = 0; i < c; ++i) {
            for (size_t j = 0; j < r; ++j){
            dividedDifferences[i*r+j] = M[1][i];
            }
        }
        for (int j = 1; j < c*r; ++j) {
            for (int i = c*r - 1; i >= j; --i) {
                if (x[i] != x[i - j]){
                    dividedDifferences[i] = (dividedDifferences[i] - dividedDifferences[i - 1]) / (x[i] - x[i - j]);
                }
                else{
                    t = t*j;
                    dividedDifferences[i] = M[j+1][i/r]/t;
                }
            }
        }

        // 构建埃尔米特插值多项式
        double result = 0.0;
        for (size_t i = 0; i < dividedDifferences.size(); ++i) {
            double term = dividedDifferences[i];
            for (size_t j = 0; j < i; ++j) {
                term *= (x1 - x[j]);
            }
            result += term;
        }
        return {result};
    }
};
class BezierInterpolation : public Interpolation {
private:
    std::vector<std::vector<double>> controlPoints; // 控制点，二维向量

public:
    BezierInterpolation(const std::vector<std::vector<double>>& controlPoints) 
        : controlPoints(controlPoints) {
        if (controlPoints.empty() || controlPoints[0].empty()) {
            throw std::runtime_error("Error: Control points must be non-empty and have at least one dimension.");
        }
    }

    double bernsteinPolynomial(int n, int i, double t) {
        double coeff = 1;
        for (int j = 0; j < i; ++j) {
            coeff *= (n - j) / static_cast<double>(j + 1);
        }
        return coeff * std::pow(t, i) * std::pow(1 - t, n - i);
    }

    std::vector<double> interpolate(double t) override {
        if (t < 0.0 || t > 1.0) {
            throw std::runtime_error("Error: t must be in the range [0, 1].");
        }

        int n = controlPoints.size() - 1;
        int dimension = controlPoints[0].size();
        std::vector<double> result(dimension, 0.0);

        for (int i = 0; i <= n; ++i) {
            double b = bernsteinPolynomial(n, i, t);
            for (int d = 0; d < dimension; ++d) {
                result[d] += b * controlPoints[i][d]; // 累加每个维度的插值
            }
        }

        return result; // 返回插值结果
    }
};
#endif // INTERPOLATION_HPP